#### FIGURE 3: NEGATIVE CONSEQUENCES (MECHANISMS)
#### Pooled effects on perceived negative consequences of immigration

rm(list = ls())
source("./2_code/00_setup.R")

#### STUDY 1 ####

data <- fread(paste0(data_path, "data_study1.csv"), header = TRUE)
data$worked <- as.factor(data$worked)

# New Conflict 
lin_conflict_welfare_economic1 <- lm_lin(new_conflict_s ~ treat_economic, covariates = ~
                                           edad + worked + male, data = data[(treat=='economic' | treat=='control')])
lin_conflict_welfare_humanitarian1 <- lm_lin(new_conflict_s ~ treat_humanitarian, covariates = ~
                                               edad + worked + male, data = data[(treat=='humanitarian' | treat=='control')])


#### STUDY 2 ####

data2 <- fread(paste0(data_path, "data_study2.csv"), header = TRUE)
data2$worked <- as.factor(data2$worked)

# New Conflict 
lin_conflict_welfare_economic2 <- lm_lin(new_conflict_s ~ treat_economic, covariates = ~
                                           edad + worked + male, data = data2[(treat=='economic' | treat=='control')])
lin_conflict_welfare_humanitarian2 <- lm_lin(new_conflict_s ~ treat_humanitarian, covariates = ~
                                               edad + worked + male, data = data2[(treat=='humanitarian' | treat=='control')])


#### STUDY 3 ####

data3 <- fread(paste0(data_path, "data_study3.csv"), header = TRUE)
data3$employed <- as.factor(data3$employed)

# New Conflict 
lin_conflict_welfare_economic3 <- lm_lin(new_conflict_s ~ treat_economic, covariates = ~
                                           edad + employed + male, data = data3[(treat=='economic' | treat=='control')])
lin_conflict_welfare_humanitarian3 <- lm_lin(new_conflict_s ~ treat_humanitarian, covariates = ~
                                               edad + employed + male, data = data3[(treat=='humanitarian' | treat=='control')])


#### META ANALYSIS ####

### Calculate observations for economic treatment ###

dataset_names <- c("data", "data2", "data3")

calculate_total_observations_eco <- function(dataset) {
  dataset %>%
    filter(treat == 'economic' | treat == 'control') %>%
    summarise(total_observations = sum(complete.cases(across(c(new_conflict_s)))))
}

total_observations <- lapply(dataset_names, function(name) {
  calculate_total_observations_eco(get(name))
})

obs_eco <- unlist(total_observations)

types_economic <- data.frame(
  treatment = rep('economic', 3),
  Study = c('Study 1', 'Study 2', 'Study 3'),
  n_total = obs_eco
)


### Calculate observations for humanitarian treatment ###

calculate_total_observations_hum <- function(dataset) {
  dataset %>%
    filter(treat == 'humanitarian' | treat == 'control') %>%
    summarise(total_observations = sum(complete.cases(across(c(new_conflict_s)))))
}

total_observations_hum <- lapply(dataset_names, function(name) {
  calculate_total_observations_hum(get(name))
})

obs_hum <- unlist(total_observations_hum)

types_humanitarian <- data.frame(
  treatment = rep('humanitarian', 3),
  Study = c('Study 1', 'Study 2', 'Study 3'),
  n_total = obs_hum
)


### Economic treatment, new conflict outcome ###

out_economic_newconflict <- rbind(
  tidy(lin_conflict_welfare_economic1) %>% filter(term=='treat_economic') %>% dplyr::select(estimate, std.error),
  tidy(lin_conflict_welfare_economic2) %>% filter(term=='treat_economic') %>% dplyr::select(estimate, std.error),
  tidy(lin_conflict_welfare_economic3) %>% filter(term=='treat_economic') %>% dplyr::select(estimate, std.error)
)

out_economic_newconflict <- cbind(out_economic_newconflict, types_economic)

m.gen_eco_newconflict <- metagen(
  TE = estimate,
  seTE = std.error,
  studlab = Study,
  data = out_economic_newconflict,
  common = FALSE,
  random = TRUE,
  method.tau = "REML",
  n.e = n_total,
  title = "economic_treatment_NEWconflict"
)

summary(m.gen_eco_newconflict)


### Humanitarian treatment, new conflict outcome ###

out_humanitarian_newconflict <- rbind(
  tidy(lin_conflict_welfare_humanitarian1) %>% filter(term=='treat_humanitarian') %>% dplyr::select(estimate, std.error),
  tidy(lin_conflict_welfare_humanitarian2) %>% filter(term=='treat_humanitarian') %>% dplyr::select(estimate, std.error),
  tidy(lin_conflict_welfare_humanitarian3) %>% filter(term=='treat_humanitarian') %>% dplyr::select(estimate, std.error)
)

out_humanitarian_newconflict <- cbind(out_humanitarian_newconflict, types_humanitarian)

m.gen_hum_newconflict <- metagen(
  TE = estimate,
  seTE = std.error,
  studlab = Study,
  data = out_humanitarian_newconflict,
  common = FALSE,
  random = TRUE,
  n.e = n_total,
  method.tau = "REML",
  title = "humanitarian_treatment_NEWconflict"
)

summary(m.gen_hum_newconflict)


#### PLOTTING ####

# Extract data for economic treatment, negative consequences outcome
total_observations <- sum(m.gen_eco_newconflict$n.e)
forest_data_eco_cons <- data.frame( 
  Study = c(m.gen_eco_newconflict$studlab, "Pooled Effect"),
  Coefficient = c(m.gen_eco_newconflict$TE, m.gen_eco_newconflict$TE.random),
  Observations = c(m.gen_eco_newconflict$n.e, total_observations),
  SE = c(m.gen_eco_newconflict$seTE, m.gen_eco_newconflict$seTE.random),
  Lower95 = c(m.gen_eco_newconflict$lower, m.gen_eco_newconflict$lower.random),
  Upper95 = c(m.gen_eco_newconflict$upper, m.gen_eco_newconflict$upper.random),
  Lower90 = c(m.gen_eco_newconflict$TE - (1.645 * m.gen_eco_newconflict$seTE), m.gen_eco_newconflict$TE.random - (1.645 * m.gen_eco_newconflict$seTE.random)),
  Upper90 = c(m.gen_eco_newconflict$TE + (1.645 * m.gen_eco_newconflict$seTE), m.gen_eco_newconflict$TE.random + (1.645 * m.gen_eco_newconflict$seTE.random)),
  Narrative = "Economic Narrative",
  Outcome = "Negative Consequences"
)

# Extract data for humanitarian treatment, negative consequences outcome
total_observations <- sum(m.gen_hum_newconflict$n.e)
forest_data_hum_cons <- data.frame( 
  Study = c(m.gen_hum_newconflict$studlab, "Pooled Effect"),
  Coefficient = c(m.gen_hum_newconflict$TE, m.gen_hum_newconflict$TE.random),
  Observations = c(m.gen_hum_newconflict$n.e, total_observations),
  SE = c(m.gen_hum_newconflict$seTE, m.gen_hum_newconflict$seTE.random),
  Lower95 = c(m.gen_hum_newconflict$lower, m.gen_hum_newconflict$lower.random),
  Upper95 = c(m.gen_hum_newconflict$upper, m.gen_hum_newconflict$upper.random),
  Lower90 = c(m.gen_hum_newconflict$TE - (1.645 * m.gen_hum_newconflict$seTE), m.gen_hum_newconflict$TE.random - (1.645 * m.gen_hum_newconflict$seTE.random)),
  Upper90 = c(m.gen_hum_newconflict$TE + (1.645 * m.gen_hum_newconflict$seTE), m.gen_hum_newconflict$TE.random + (1.645 * m.gen_hum_newconflict$seTE.random)),
  Narrative = "Humanitarian Narrative",
  Outcome = "Negative Consequences"
)


### Combined plot ###

forest_data_combined_cons <- rbind(forest_data_eco_cons, forest_data_hum_cons)

desired_order <- c("Pooled Effect", "Study 3", "Study 2", "Study 1")

forest_data_combined_cons <- forest_data_combined_cons %>%
  mutate(Study = factor(Study, levels = desired_order))

coefficient_plots_cons <- ggplot(forest_data_combined_cons, aes(x = Coefficient, y = Study)) +
  geom_point(aes(color = Study), size = 3) +
  geom_errorbarh(aes(xmin = Lower95, xmax = Upper95, color = Study), height = 0, linewidth = 0.9) +
  geom_errorbarh(aes(xmin = Lower90, xmax = Upper90, color = Study), height = 0, linewidth = 1.9) +
  geom_vline(xintercept = 0, linetype = "dashed", color = "red") +
  labs(title = "", x = "", y = "", caption = "ATE Negative Consequences") +
  theme_minimal() +
  theme(
    axis.title.x = element_blank(), 
    plot.caption = element_text(size = 15, hjust = 0.5),
    strip.text = element_text(size = 17),
    axis.text.y = element_text(size = 13),
    axis.text.x = element_text(size = 13),
    legend.position = "none"
  ) +
  scale_x_continuous(limits = c(-0.3, 0.1)) +
  scale_color_manual(values = c("Study 1" = "grey50", "Study 2" = "grey50", "Study 3" = "grey50", "Pooled Effect" = "red")) +
  facet_grid(~Narrative)

coefficient_plots_cons

ggsave(filename = paste0(plot_path, "figure_03.png"), plot = coefficient_plots_cons, width = 10, height = 4.5, dpi = 600)